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Abstract 

A model of three-body motion is developed which includes the effects 
of gravitational radiation reaction. The radiation reaction due to the 
emission of gravitational waves is the only post-Newtonian effect that is 
included here. For simplicity, all of the motion is taken to be planar. 
Two of the masses are viewed as a binary system and the third mass, 
whose motion will be a fixed orbit around the center-of-mass of the bi- 
nary system, is viewed as a perturbation. This model aims to describe the 
motion of a relativistic binary pulsar that is perturbed by a third mass. 
Numerical integration of this simplified model reveals that given the right 
initial conditions and parameters one can see resonances. These (m,n) 
resonances are defined by the resonance condition, mui = 2nf2, where m 
and n are relatively prime integers and u; and are the angular frequen- 
cies of the binary orbit and third mass orbit (around the center-of-mass 
of the binary), respectively. The resonance condition consequently fixes 
a value for the semimajor axis of the binary orbit for the duration of the 
resonance; therefore, the binary energy remains constant on the average 
while its angular momentum changes during the resonance. 
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1 INTRODUCTION 

In 1975, Hulse and Taylor found through the timing observations of binary pul- 
sar PSR B1913+16 that the semimajor axis decayed at nearly the rate predicted 
by general relativity for the emission of gravitational radiation. This has served 
as significant indirect evidence of gravitational radiation reaction force as the 
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result of energy balance due to the fact that gravitational waves leave the sys- 
tem. It has also shown the importance of the pulsar in probing strong gravity. 
Because of the universal nature of gravity, a correct analysis of astrophysical 
systems must involve the gravitational interaction between all of the masses in 
the environment plus other non-Newtonian effects. One would expect that over 
time, evidence of neighboring masses would appear in the timing observations 
of binary pulsars. 

This study puts forth a model of three bodies, two of which constitute a rel- 
ativistic binary pulsar and the third is a perturbing mass. This investigation ex- 
tends the work of Chicone, Mashhoon, and Retzloff (1996a,1996b,1997a,1997b), 
which originally dealt with a binary system perturbed by normally incident 
gravitational waves. As in these previous models, the model derived here will 
be searched for resonances. In the following analysis, an (m, n) resonance will 
occur when the relation mu> = 2nfl is satisfied, where m and n are relatively 
prime integers, co is the frequency of the binary system, and Q is the frequency 
of the orbit of the third mass around the center-of-mass of the binary. Such 
resonances are physically noteworthy because they correspond to events where 
the collapse of the semimajor axis of the binary halts. That is, on average it 
stays fixed. 

The N-body problem in general relativity theory is rather complicated; there- 
fore, to capture the main effects of three bodies plus gravitational radiation 
reaction one can start with a classical three-body system with gravitational 
radiation damping as a linear perturbation. That is, in this work all post- 
Newtonian (relativistic) corrections will be neglected except for gravitational 
radiation reaction that will be treated in the quadrupole approximation. For 
the sake of simplicity, the internal structure of the masses will be neglected so 
that one is in effect dealing with three Newtonian point masses. 

Two astrophysical systems that may demonstrate measurable behavior based 
on this model are the relativistic binary pulsar near a massive third mass and 
a planet in an orbit around a binary pulsar. While this study will focus on the 
former case the latter proves to be a timely point of interest. Astronomers have 
in the last six years found more than 50 stars which include at least one planet. 
The present number of stars that are currently being observed for planets is 1000 
(Cameron, 2001). Pulsars have also proven to be detectable sources that show 
evidence of extrasolar planets. Finding a pulsar planet, where the pulsar is also 
part of a binary system is a possibility. Microlensing technique in its extrasolar 
planet search has turned up a possible planet with two suns (Bonnet et al., 
1999). Another possibility would be the future detection of a binary pulsar in 
a globular cluster. A binary pulsar system that is situated in an approximately 
spherical cluster would be effectively the same as a binary attracted to a third 
mass equal to the mass of the matter inside the binary's orbit. 

The timing of radio pulses from pulsars has opened up a rich and exciting 
branch of astronomy and astrophysics. Analysis of pulsar signals over time offers 
insight into the gravitational environment such as the possible presence of a 
companion mass. Hence, pulsar astronomy offers a fruitful method of detection 
for the astrophysical processes discussed in this paper. On the frontiers of 
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astronomy, interferometric gravitational wave detectors based on Earth such 
as LIGO, VIRGO, GEO, AND TAMA aim to detect gravitational waves from 
known inspiraling binary systems. The shortening orbital separation that would 
arise as the result of energy loss via gravitational wave emission would result 
in an increasing frequency in the orbit. The gravitational wave signal would 
have twice this frequency (Blanchet, 2002). Thus, the predicted effect of a 
binary system's capture into resonance would leave its imprint in terms of the 
detected gravitational waves. Namely, the interval when the distance between 
the members of the binary stays fixed on average would give a gravitational wave 
signal that on the average would be of near constant frequency. Furthermore, 
the next generation of planned space-based detectors like LISA will broaden the 
scope of gravitational wave detection. 

The initial approach to this problem in Section 2 involves setting up the 
equations of motion for the three masses with the influence of gravitational 
radiation reaction included. The radiation reaction force is understood to be 
very small, so it is viewed as a perturbation. The model is then simplified in 
Section 3 to make it more amenable to analysis. The results of the numerical 
integration of the simplified model are presented in Section 4. Next, section 5 
expands upon the numerical results with a discussion about resonance. Section 
6 contains concluding remarks. 

The nonlinear disposition of the equations of motion for the relative motion 
of the binary system requires analysis that is suited for nonlinear behavior. 
Chicone, Mashhoon, and Retzloff (1997b) developed such an averaging method, 
that elucidates details about the structure of the orbit of a nonlinear system , 
especially when a resonance occurs. The system developed and studied here is 
suited for further analysis, in particular the application of the method mentioned 
above. 



2 EQUATIONS OF MOTION 

The general equations of motion for the three-body problem with gravitational 
radiation damping are: 

d 2 x\ Gm 1 m 2 (x\ - x\) _ Gmim 3 (x\ - x l 3 ) 2G d 5 D lj j 
mi ^ + | Xl -x 2 |3 = | Xl -x 3 |3 ' "lS?™ 1 "^-* 1 ' W 

d 2 x% Gm\m 2 (x\ — x\ ) Gm 2 m 3 (x% — x\) 2G d 5 Da ,■ 

m2 ^7T + 1 i3 = 1 | 3 - T^rn 2 — -^x> 2 , 2 

dt z l x i — x 2| |x 2 — X3| d 15c 5 dt b 

d?x 3 Gm 3 m 1 (x 3 — x\) Gm 2 m 3 (x 3 — x\) 2G d 5 Dij j 

m3 ^ + |X3-X!|3 = |X 3 -X 2 |3 l^™ 3 ^^ 3 ' (3) 

where the quadrupole moment tensor for the three-body system is: 

Dij = mi(3x\x{ — 8ijx\) + m2(ix l 2 x 3 2 — 5ijX 2 ) + 7^3(3x3X3 — SijX s ). (4) 
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Each equation includes the familiar Newtonian gravitational interactions 
between the masses plus the Newtonian gravitational radiation reaction force. 
This damping force results from the emission of gravitational waves by the 
system of masses. There is a loss of energy in the system that results from 
energy carried away by the waves. According to general relativity, the energy of 
gravitational waves leaving the system is given in the quadrupole approximation 
by 

dE _ G d 3 D t3 d 3 D^ 

~dt ~ 45^~dfi dW' ( ' 

where E is the amount of gravitational radiation energy emitted by the system 
(Landau and Lifshitz, 1971). 

Energy conservation requires that this loss be reflected in the equations of 
motion; therefore, one can also arrive at this rate of energy loss by way of 
mechanics. If one multiplies each equation of motion by its time derivative of 
position and adds the three equations, one arrives at the following expression: 

df L _ G d 3 D i:j d 3 D^ 

~dt ~ _ 45?~^3 dt^ 1 ( ' 

where £ = £n + £s- Here £n is the total Newtonian energy of the three-body 
system. That is, 

i=i ijtj J 

and £s is the gravitational "Schott" energy, 

G ( dtDijdPv d-i),, ,/ '/)'■•- \ 
S 45c 5 ^ dt 4 dt dt 3 dt 2 J ' [ ' 

that is analogous to the Schott term in electrodynamics (Schott, 1912). 

For quasi-periodic motions, the change in the Schott energy over a period 
turns out to be of higher order and can be neglected. Hence there is agreement 
between the two different methods as to what the loss of energy should be 
since the right-hand side of (5) is equal and opposite to the right-hand side of 
(6). Comparison of the general relativistic and classical methods of calculating 
rate of change of angular momentum yields agreement in a similar way. Thus 
in the quadrupole approximation for the emission of gravitational waves, the 
three-body system loses energy and angular momentum to the radiation field. 
Moreover, adding equations (l)-(3) and making the following definition 

z _ mixi + m 2 x 2 + to 3 x 3 
mi + m 2 + m 3 

for the center-of-mass of the system we find 

d 2 Z i 2G d 5 Dij ■ 

^r + i5^^r z =0 - (10) 



4 



Let us note that Z = is a solution of this equation, so that if the center-of-mass 
of the whole system is initially at rest at the origin of coordinates, it will remain 
so in the quadrupole approximation under consideration here. This is consistent 
with the fact that gravitational waves do not carry away linear momentum in 
the quadrupole approximation. 

Solutions to the proposed three-body problem, equations (l)-(3), by way of 
numerical integration offer some insight into the behavior of the system that 
has no known solution in closed form. To integrate such a higher-order system 
numerically, it is not sufficient to specify the positions and velocities of the 
three masses at some initial instant of time. In fact, it is known that such 
systems suffer from the existence of runaway modes that inevitably will lead 
to divergent results. It is therefore necessary to replace the system (l)-(3) by 
second-order equations of motion via iterative reduction as explained in a recent 
paper by Chicone, Kopcikin, Mashhoon, and Retzloff (2001). This procedure 
involves the repeated substitution of the equations of motion (l)-(3) in the 
evaluation of d 5 Dij/dt 5 and the subsequent reduction of the resulting system 
to one of second-order equations for the motion of the system (by dropping 
higher-order terms). The resulting system would be appropriate for numerical 
integration; however, it would have a rather complicated form. Therefore, we 
resort to certain simplifications in this first treatment of the classical three-body 
problem that takes gravitational radiation reaction into account. 

The main physical model which will be under investigation in this study 
involves two masses, which represent the binary, whose center-of-mass orbits a 
massive third body. Since the mass of the third body is assumed to be much 
larger than the mass of the binary, one can take X3 to be the origin of the 
coordinate system. This is because the center-of-mass of the entire system is 
very close to the location of the third mass. To simplify matters, we set x 3 = 0; 
moreover, the distance of the binary from the third mass is taken to be much 
larger than the semimajor axis of the binary system. 

This model could approximate a relativistic binary pulsar system that orbits 
the center-of-mass of a globular cluster. The outer shell of stars would remain 
'unseen' by the binary since a globular cluster is nearly spherical, and the interior 
stars would serve as a third mass located at the center-of-mass of the cluster. 



3 SIMPLIFIED MODEL 

The following equations of motion result from the approximations made in the 
model: 

d 2 x\ Gm 3 mix\ Gm\n\2{x\ — x\) 2rriiGd 5 Dij j 
mi ^ = \^\^ + |x 2 - Xl |3 15J~~W~ Xl ' ( > 

d 2 x\ Gm 3 m2X l 2 Gm,\m,2{x\ — x\) 2ni2Gd 5 Dij j 

™ 2 ^~ = ^ |x 2 -Xi|3 15^^" ^ 

Here mi and m<2 are masses of the members of the binary, 7713 is the mass of 
the perturbing third body, and Dij = mi(3x\x{ — Sijxf) + m2(3x2X 2 — SijX 2 .) 
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is the quadrupolc moment tensor. 

In order to elucidate the relevant dynamics of the system, it is helpful to 
change the coordinates of the system. The relative and center-of-mass coordi- 
nates of the binary can be defined as: 

rW 2 -*i, * = "*'i (13) 

where M = mi + m 2 . The equation of relative motion is then 

rfV _ GMr* /_4 x\ \ 2G d 5 D tJ ■ 

dt 2 |r| 3 3 ^|x 2 | 3 | Xl |3j 15c 5 dts r ' 

while the equation of motion of the center-of-mass of the binary is 

d 2 X % Gm?,(m\x\ m 2 x 2 \ 2G d 5 D i:j 

' X J . (15) 



A 2 M V Ixil 3 |x 2 |3 J 15c 5 A 5 



One can show that the quadrupolc moment tensor can be written in terms of 
the relative and center-of-mass variables as = fj,(3r l r 3 — Sijr 2 ) + M(3X 1 X : > — 
5ijX 2 ), where fi = mim 2 /M is the reduced mass of the binary system, r = |r|, 
and X = |X|. Replacing xi and x 2 in (14)-(15) with x\ — X 1 — 72 r l and 
x\ = X 1 + "fir 1 and expanding terms in powers ofr/X, one obtains 



X 1 - 72 r* X 1 72 / . 3(X • r)X 



|xi| 3 |X- 72 r| 3 X 3 X 3 V X 2 
xi X l + 7l r* X 1 7l / . 3(X-r)X< 



|x 2 | 3 |X + 7 ir| 3 X 3 X 3 V X 2 



+ ..., (16) 
+ ••., (17) 



where 7 i = m\/M and 72 = m 2 /M. 

The following equations result when terms to linear order in r/X are kept 
in the above expansions: 

dV _ _GM Grr^ 2G d^j , 

dt 2 ~ r 3 X 3 13 15c 5 dtB ' U8j 

W ^.i^d^ (19) 



where 



(it 2 X 3 15c 5 dt 5 



is the reduced tidal matrix. 

Thus the problem reduces to the coupled system (18)-(19) for the motion 
of the binary system. To simplify the equations further, it will be assumed in 
what follows that the center-of-mass motion is circular in the absence of radi- 
ation reaction. Substituting such a solution for (19) in the equation of relative 
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motion (18), one finds that to lowest order in the perturbations the center-of- 
mass motion can be taken to be simply circular in (18). This circular motion 
contributes to the radiation reaction term in (18); however, this contribution 
can also be neglected if 

X\ n / 2 HH/ m 1 + m 1 \ 1/2 >:> L 
r J \m 3 J \m 3 J \ m 3 J 

Assuming this inequality, the equation of relative motion reduces to 
dV GM • Gm 3 „ . 2Gfid 5 D l 



-r — 



K^--^—^r\ (22) 



dt 2 r 3 X 3 15c 5 dt 5 

where Dij = 2>r % r^ — Sijr 2 . The rest of this paper is devoted to the study of this 
equation that describes the relative motion of two bodies under mutual gravi- 
tational attraction, tidal interaction with a large third mass, and gravitational 
radiation damping. These are the combined effects of interest here. This model 
can, in effect, represent a second scenario as well, i.e. a binary system orbited 
by a distant small mass m 3 (see Appendix A). 

It is convenient to transform equation (22) into dimensionless form. To 
this end, let all lengths and temporal variables be measured in units of Ro 
and To, respectively. Moreover, one assumes that these units are related to 
the unperturbed motion of the binary such that GMTq = Rq. Thus, letting 
r % — > R f l , X 1 — » RoX 1 , and t — ► T i in equation (22), it reduces to the form 

d2fl _ fl cfK r 1 5 d5Dij r> (2^) 
dP ~ f 3 13 dP ' 1 j 

where the hats will be dropped in what follows for the sake of simplicity. Here 

3 

Y • ~ A I \ X 
2G^i?| 

d ~15c5T3' [2b > 

and the center-of-mass motion is taken to be a circle in the (x, y)-plane such 
that (Kij) has the form 



V = ^ = ^(^) , (24) 



-i-|cos20i -| sin 2Qt 
-| sin 2nt -| + | cos 2flt 

1 



(26) 



If the length unit i?o were chosen to be the semimajor axis of the binary, 
the corresponding initial period would be 27rT . A relativistic binary pulsar 
such as PSR B1913+16 has an orbital period of nearly eight hours (Lyne and 
Graham-Smith, 1998). The data regarding this Hulsc- Taylor binary pulsar when 
applied to the dimensionless formula (25) give 5 ~ 1CT 16 . This value of S is also 
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approximately valid for the relativistic binary pulsar PSR B1534+12 (Stairs et 
al., 1998). Moreover, the tidal perturbation term is also assumed to be very 
small, i.e. ft 2 « 1. It proves useful to replace f2 2 by a free parameter e << 1 
in what follows. Therefore, equation (22) will be replaced by 

dt 2 r 3 1 dt 5 K ' 

The mathematical results derived from (27) will then apply to the physical 
situation at hand once e = ft 2 . 

Finally, the order of equation (27) needs to be reduced in order to avoid a 
singular perturbation problem. To this end, in taking derivatives of Dij substi- 
tutions for P are made with the equation of motion (27), and terms to desired 
order are kept. In this case, we reduce (27) to a second-order equation that is 
linear in e and 5. See Appendix B for a complete derivation. What results is 

— = --- eKi ^-5R\ (28) 

Here W is the reduced radiation reaction term in Cartesian coordinates and is 
given by 

Ri = (-y - 18W2 + 72w 2 ) ^ + (-* + 300f 2 - 216« 2 ) ^J, (29) 
where as before r = Irl, v % — f l and v = Ivl. 



4 NUMERICAL RESULTS 

What results is a nonlinear two-dimensional second-order ODE that describes 
the relative motion between the two members of the binary. Because of the non- 
linearity, a numerical solution is sought at this stage to elucidate the behavior 
of the system. The numerical integrator used in this study was furnished by 
MATHEMATICA. To produce a trajectory of the system, one needs values for 
the parameters e, 5, and plus the necessary initial conditions. For the pur- 
poses of generating numerical solutions for inspection, it becomes advantageous 
to generalize the system as a dynamical system at the expense of adherence to 
a strictly physical model. This relaxation of the parameter space allows one to 
more freely search for resonances in the generalized system. 

Because near-Keplerian motion is under investigation, a natural analytical 
convention to use is the osculating ellipse. The otherwise constant elements of 
the ellipse for unperturbed motion become time-dependent osculating elements 
with the perturbations. Each instantaneous position in the motion is a point of 
an ellipse that is characterized by the osculating elements at that time. If at that 
time the perturbation were supressed, the motion would follow the path of the 
osculating ellipse uniquely described by the instantaneous position and velocity. 
One can, of course, transform between coordinate systems. The scheme for 
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numerical work, for example, was done in polar coordinates. When one takes 
equation (28), applies the transformation to polar coordinates 

x = r sin 0, y — r cos 9, (30) 

and converts the second-order ODEs into equivalent first-order ODEs one ob- 
tains the following system of equations: 

r = P r , (31) 
8=£, (32) 

p2 i i p oo p2 

P r = -§--- + -er[l + 3cos (n't - 29)] + <5-f (- + 24P 2 + 144-f-), (33) 
P e = -er 2 sin (fi't - 29) +*-£(- + 108P 2 - 72-J-). (34) 

Z (I 1 

Here f2' = 20; furthermore, this polar scheme with the usual variables r and 6> 
gives rise to the variables P r and Pe defined above in (31) and (32) in this first 
order form. Pq represents the angular momentum of the two-body system per 
reduced mass. 

One can generate a series of numerical integrations in search of resonances of 
various orders. Once the parameters that represent the amplitudes of the tidal 
perturbation due to the presence of the third mass and gravitational radiation 
damping arc chosen, one can numerically generate an orbit starting from initial 
conditions. In anticipation of future analytical work to be done in Delaunay 
variables (Wardell, to be published), results of the numerical integration are 
graphed in terms of Delaunay variables. In figure 1, the orbit of the Delaunay 
action L is plotted versus time. The variable L is directly related to the semima- 
jor axis by the relation L = a 1 / 2 , where a is the semimajor axis of the relative 
orbit. The left side of the orbit makes a sharp descent which accords with the 
expectation of semimajor axis decay due to the influence of gravitational ra- 
diation reaction. However, the orbit commences to undergo oscillations about 
an average value. A resonance occurs here. The amplitude of the oscillations 
increase until the orbit falls out of resonance. This is a nonlinear feature of the 
orbit. The orbit eventually falls out of the resonance and gives way to a descent 
whose slope is more gradual than that of the initial sharp descent due to the 
fact that the orbital eccentricity changes during resonance as a consequence of 
the variation of orbital angular momentum. The phenomena associated with 
this resonance are discussed in the next section. 

5 RESONANCE 

Numerical experiments show that a (1 : 1) resonance exists for the system with 
the appropriate initial conditions. Checks on the numerical integration have 
been performed to rule out the possibilty that the effects seen are numerical 
artifacts. 
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The gravitational environment of the binary comes from two contributors: 
the tidal influence of the third mass and the emission of gravitational waves. 
These combined effects enter as the perturbations in the dynamical system. The 
interplay of tidal energy input from the orbit around the third mass and the 
energy loss caused by the emission of gravitational waves comes to a place of 
balance during a resonance. A resonance is characterized by the resonance con- 
dition mu> = 2nf2, where u> is the angular frequency of the perturbed Kepler 
orbit and fl is the angular frequency of the binary's orbit around the third mass. 
Because the radiation damping is a dissipative effect, decay of the semimajor 
axis might be expected as the dominant pcrturbativc effect. In fact, the exam- 
ple given by the Hulse- Taylor binary pulsar where this decay is predominantly 
observed shows this effect quite well. However, in the presence of a third mass 
resonance, as seen in this system, fixes the average net flux of energy and gives 
rise to an interesting nonlinear dynamical system whose perturbations of tidal 
interaction and gravitational radiation damping offset each other on the aver- 
age. Changes in the orbital angular momentum, though, can accompany this 
fixed average energy that follows from the resonance condition. The binary's 
orbit around the third mass can result in an increase in its internal orbital an- 
gular momentum due to the tidal torque of the external mass. On the other 
hand, gravitational waves emitted from the binary carry with them angular 
momentum that decreases the orbital angular momentum in the binary. 

5.1 CAPTURE INTO RESONANCE 

Approaching from the left in Figure 1, the semimajor axis of the orbit steadily 
shrinks until it reaches the resonance manifold. The resonance manifold is de- 

— 1/3 

fined by L = L„, where L* = ui* and the resonance condition for a (m : n) 
resonance fixes such that uj* = 2(n/m)£l. This in turn fixes L*. Resonance 
capture is noteworthy because the orbit of the system passes through many reso- 
nances unimpeded. The integers m and n must be relatively prime, making n/m 
a rational number; since the rational numbers are dense with respect to the real 
numbers, the orbit is always close to a resonance without being captured, except 
under special circumstances. Numerical experiments show that capture more 
readily occurs for low- lying resonances such as (1 : 1) as illustrated in Figures 
1-3. Moreover, higher order resonances give rise to more chaotic structures in 
their corresponding figures for L, Pg, and e (Chicone, Mashhoon, and Retzloff, 
1997b). Here Pg is the specific orbital angular momentum of the osculating 
ellipse and e is its eccentricity defined by (1 — Pg /L 2 ) 1 / 2 . 

5.2 PASSAGE THROUGH RESONANCE 

When an orbit is captured into a resonance it enters the resonance manifold. 
That is, the average L value is fixed at L*. The orbit then oscillates about this 
value with increasing amplitude. Eventually, it falls out of the resonance. 
While in resonance, the other action variable Pg can on average change signifi- 
cantly. This is shown in Figure 2, which indicates that through the resonance L 
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oscillates about an average value L», Pg increases on average during this inter- 
val. Thus in the course of resonance, the osculating ellipse associated with the 
orbit undergoes an average change in eccentricity, as demonstrated in Figure 3. 

5.3 EXIT FROM RESONANCE 

The energy loss due to the emission of gravitational waves depends on the 
eccentricity as (1 — e 2 )~ 7//2 (Landau and Lifshitz, 1971). Therefore, the decrease 
in this orbital eccentricity accounts for the apparent lower rate of energy loss 
once the orbit leaves the resonance as in Figure 1. The anti-damping seen 
through the resonance eventually disrupts the resonance condition and leads 
to the orbit falling out of resonance. Where the orbit falls out of resonance 
the eccentricity is e ~ 0.73 as compared to e ~ 0.85 when the orbit was first 
captured into resonance. 

An approximate analytic description of the phenomenon associated with 
resonance is the subject of a future publication (Wardell, to be published). 

6 CONCLUSION 

The main thrust of the analysis in this paper has been to introduce and develop 
a model of a binary system that is under the influence of gravitational radiation 
reaction and perturbed by a third body. The classical case of three gravitating 
bodies is considered here when the damping force due to the emission of gravi- 
tational waves by the system is included in the Newtonian equations of motion. 
After the pertinent equations of motion were derived, numerical analysis was 
done to explore important dynamics in the system - namely, resonance between 
the orbits of the relative and third body motions. This would optimistically be 
an astronomically observable effect. 

An analytical approach that would resolve details about this nonlinear sys- 
tem would add to the confidence of the numerical results. What remains, to 
augment the results of this work, is to pursue an analytical solution. The pre- 
vious work of Chicone, Mashhoon, and Retzloff (1996,1997) outlines a novel 
averaging approach that looks into orbits near resonance. It was applied to the 
case of a binary system perturbed by the emission and absorption of gravita- 
tional waves. The model presented in this paper along with numerical indication 
of the existence of resonances leads naturally to this analytical approach for the 
model developed here (Wardell, to be published). 
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Figure 1: This is a (1:1) resonance. The parameters are e = 0.00003, 5/e = 
10~ 3 , and Q, = 0.14286. The initial conditions at t = are (r,0,P r ,P e ) = 
(1,0.57,0.75246,1). 
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Figure 2: Plot of Pg versus time for the (1:1) resonance. 
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Figure 3: Plot of eccentricity versus time for the (1:1) resonance. 
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APPENDIX A: THE SECOND SCENARIO 

Imagine a relativistic binary system consisting of masses mi and and a 
distant third mass m 3 (with to 3 << m\, and m 3 << m-i) in a nearly circular 
orbit about the binary system such that r << |x 3 |. Starting from (l)-(4), the 
equation of relative motion is 



dt 2 + 3 ^|x 3 -x 2 |3 |x3-xi|V 15c 5 dt* [ ' 

Writing x\ = X 1 — and x\ = X 1 + 71 r % as before, and recalling that 
X 1 = — m^x\/M, we find 



(l + ^)4- 7 ir' 

|3 in J- 223. Wn - -v,r|3 : 



(A2) 



|X3-X 2 |3 |(l + ^)x 3 -7ir| 
)xl 4- ^r' 

(A3) 



4-*j = (l + ^)4+72r' 

|X 3 - X1 |3 |(l + IM )x3+72r |3- 

Using m 3 << M and neglecting 777,3 /M in comparison to unity meanwhile ex- 
panding (A2) and (A3) as in equations (18) and (19) based on the fact that 
r << |xa|, one finds that to lowest order the result is 

dV GM _ Grn, 2G_ d^ 



dt 2 r 3 |x 3 | 3 J 15c 5 dt 5 
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where ~ ^(3rV J — %r 2 ) once the condition 

/ \ n/ 2 

M VI X 3|/ 

is imposed. In equation (A4), Kij is given by 

Kii=5ij-Z0- (A6) 

It is clear that this scenario, suitably interpreted, gives essentially the same 
equation of relative motion that is given in equation (27). 

APPENDIX B: ITERATIVE REDUCTION 

To obtain the desired approximation of the radiation reaction force, one makes 
successive differentiations of the quadrupole moment tensor and substitutes the 
Keplerian equation of motion whenever the second time derivative of position 
appears. The reduced radiation reaction is: 

W - (Bl) 



The scaled Kepler equation of motion is: 



r 

(B2) 



dt 2 r 3 ' 

For brevity, r will be represented by v. One starts with the quadrupole moment 
tensor 

D i:j = Sr'r 3 - 8 l3 r 2 (B3) 
and proceeds to take time derivatives: 



dt 

d 2 D^ 



<lD ii = 3 ( r v + «V) - <% (2r • v) , (B4) 



f- = 3 (rV + t)V + 2vV) - 2S i3 (v 2 + r • v) . (B5) 



dt 2 

Where v l occurs, one substitutes the unperturbed Kepler equation of motion to 
arrive at the following 

*>-«(-£+«v)- ^-i). (B» 



Subsequent differentiations and substitutions yield the following expressions: 
^ = ~ («V + rV) + 18^r - % (-2^) , (B7) 
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V = " — + ^3 ^ + (-90^ + 18« 2 + - ) — (B8) 

S f^-r 2 -— — 

d 5 A, 18(W' . / 24 „„ ,\ /VW + rV 

rfi 5 r 



+ (630f 2 - 270v 2 ) — r - 5- 



r + 360r^ + 72?/ = (B9) 



rV . „ / 30r 3 I8v 2 r 16r 



r 6 u ^ \ ^,4 ,,4 r 5 



Finally, when one contracts expression (B9) with to compute the desired 
reduced radiation reaction expression as prescribed by (Bl), one finds equation 
(29). 
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